Wstęp

W niniejszym raporcie przedstawiono działania podjęte w celu określenia głównych przyczyn zmniejszania się długości śledzi oceanicznych wyławianych w Europie. Predykcja została dokonana w oparciu o model przygotowany na zbiorze danych opisanym w sekcji Zbiór danych. W oparciu o przyjęte miary najlepszy okazał się model z algorytmem Random Forest.

W wyniku przeprowadzonej analizy jako czynniki mające wpływ na zmniejszenie się długości śledzi wskazano zmianę temperatury przy powierzchni wody, a także dostępność planktonu Calanus finmarchicus gat. 1 oraz widłonogów gat. 1.

Ustawienia środowiska

Biblioteki

Do przeprowadzenia analizy oraz wygenerowania raporu zostały wykorzystane następujące biblioteki:

library('knitr')
library('dplyr')
library('corrplot')
library('ggplot2')
library('randomForest')
library('caret')

Zapewnienie powtarzalności

Celem zapewnienia powtarzalności oraz odtwarzalności zaprezentowanych wyników analizy ustawiono stałą wartość ziarna generatora (ang. seed). Fragment kodu odpowiedzialny za tę funkcjonalność został przedstawiony poniżej.

set.seed(123)

Zbiór danych

Dane wykorzystane do analizy stanowią pomiary śledzi oraz warunków w jakich żyją z ostatnich 60 lat. Dane te zostały zebrane podczas połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Zbiór danych został pobrany ze strony http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/sledzie.csv (dostęp 10-11-2016).

Wczytanie danych

W analizowanym zbiorze danych znajdują się wartości puste, nieznane (NA, ang. Not Avaible), które zostały w tym przypadku oznaczone symbolem ?.

herr_data <-read.csv("sledzie.csv", na.strings = "?", comment.char = "")

Opis atrybutów

Znaczenie kolejnych kolumn reprezentujących atrybuty zostały przedstawione w poniższej tabeli.

atrybut informacja jednostka
length długość złowionego śledzia cm
cfin1 dostępność planktonu Calanus finmarchicus gat. 1 zagęszczenie
cfin2 dostępność planktonu Calanus finmarchicus gat. 2 zagęszczenie
chel1 dostępność planktonu Calanus helgolandicus gat. 1 zagęszczenie
chel2 dostępność planktonu Calanus helgolandicus gat. 2 zagęszczenie
lcop1 dostępność planktonu widłonogów gat. 1 zagęszczenie
lcop2 dostępność planktonu widłonogów gat. 2 zagęszczenie
fbar natężenie połowów w regionie ułamek pozostawionego narybku
recr roczny narybek liczba śledzi
cumf łączne roczne natężenie połowów w regionie ułamek pozostawionego narybku
totaln łączna liczba ryb złowionych w ramach połowu liczba śledzi
sst temperatura przy powierzchni wody °C
sal poziom zasolenia wody Knudsen ppt
xmonth miesiąc połowu numer miesiąca
nao oscylacja północnoatlantycka mb

Podstawowe informacje

Dane są zorganizowane wierszowo - każdy wiersz odpowiada pojedynczemu połowowi, natomiast kolumny odpowiadają zarejestrowanym zmiennym (atrybutom).

Zbiór danych składa się z 16 atrybutów oraz 52582 rekordów.

kable(summary(herr_data))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000 Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849 Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778 1st Qu.: 2.469 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808 1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:13.60 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
Median :26291 Median :25.5 Median : 0.1111 Median : 0.7012 Median : 5.750 Median :21.673 Median : 7.0000 Median :24.859 Median :0.3320 Median : 421391 Median :0.23191 Median : 539558 Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
Mean :26291 Mean :25.3 Mean : 0.4458 Mean : 2.0248 Mean :10.006 Mean :21.221 Mean : 12.8108 Mean :28.419 Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973 Mean :13.87 Mean :35.51 Mean : 7.258 Mean :-0.09236
3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936 3rd Qu.:11.500 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232 3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:14.16 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958 Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736 Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000
NA NA NA’s :1581 NA’s :1536 NA’s :1555 NA’s :1556 NA’s :1653 NA’s :1591 NA NA NA NA NA’s :1584 NA NA NA

Przygotowanie zbioru danych

Brakujące dane

Kolejny etap analizy stanowi oczyszczenie danych, którego istotnym elementem jest przetworzenie brakujących danych. Po wnikliwej analizie surowego zbioru danych zauważono, że brakujące wartości można uzupełnić na podstawie grup utworzonych na podstawie wartości atrybutów totaln oraz recr. W ramach tak zdefiniowanych grup wartości posczególnych atrybutów z brakującymi wartościami są takie same. Aby dodatkowo uniezależnić się od wpływu ewentualnych wartości odstających w ramach grupy (ang. outliers) brakujące wartości są estymowane za pomocą mediany.

na_cols <- names(which(colSums(is.na(herr_data))>0))

fill_na <- function(x) ifelse(is.na(x), median(x, na.rm=TRUE), x)
herr_data_nona <-  herr_data %>% 
  group_by(totaln, recr) %>% 
  mutate_each(funs(fill_na), one_of(na_cols)) %>%
  ungroup

Analiza korelacji pomiędz atrybutami

Następną fazę przetwarzania stanowiła analiza korelacji pomiędzy zmiennymi. Wartość współczynników korealcji pomiędzy poszczególnymi atrybutmi został zwizualizowany na poniższym rysunku.

herr_mcorr <- cor(herr_data_nona)
corrplot(herr_mcorr, type = "upper", method="color", addCoef.col = "black", tl.col="black", number.digits=2, number.cex=0.75)

Z analizy powyższej macierzy korelacji wynika, że zachodzi duża korelacja między zmiennymi chel1 oraz lcop1 (0.95), chel2 oraz lcop2 (0.89), a także fbar oraz cumf (0.82). Wobec występowania wysokiej korelacji pomiędzy wspomnianymi zmiennymi atrybuty lcop1, chel2 oraz cumf zostały wykluczone, jako redundante, nie wnoszące nowych informacji do analizy.

herr_data_nona <- herr_data_nona[ , -which(names(herr_data_nona) %in% c("lcop1", "lcop2", "cumf", "xmonth"))]

Histogram

for(i in 2:5) {
  print(
    ggplot(herr_data_nona, aes(x = herr_data_nona[i])) + 
      geom_histogram(bins = 30, fill="#56B4E9") + 
      labs(x = colnames(herr_data_nona)[i], y = "liczba wartości") + 
      ggtitle(paste("Histogram atrybutu", colnames(herr_data_nona)[i])) + 
      theme_light()
  )
}

Zmiana długości śledzi w czasie

Poniższy wykres przedstawia zmianę długości śledzi w czasie (zgodnie z opisem dane były uporządkowane chronologicznie). Niestety w samym zbiorze, jak i jego opisie nie zostały wyspecyfikowane dokładne ramy czasowe pomiarów.

library("plotly")

plot <- ggplot(herr_data_nona, aes(x=X, y=length)) + 
  geom_point() + 
  geom_smooth(method="auto", se=TRUE, color="red") +
  labs(title="Wykres zmiany długości śliedzi w czasie", x="Czas", y="Długość śledzia")

ggplotly( plot )

Wykres potwierdza postawiony problem badawczy malejącej od pewnego momentu długości śledzia.

Regresor

Podział zbioru danych

W celu stworzenia i oceny regresora konieczne jest podzielenie zbioru danych na dane uczące, walidujące i testowe. Dane uczące i walidacyjne zostały wykorzystane do budowy i optymalizacji parametrów modelu. Natomiast dane testowe posłużyły do wyboru modelu oraz jego końcowej oceny.

in_training <- 
    createDataPartition(
        # atrybut do stratyfikacji
        y = herr_data_nona$length,
        # procent w zbiorze uczącym
        p = .75,
        # zwracamy indeksy
        list = FALSE)

training <- herr_data_nona[ in_training,]
testing  <- herr_data_nona[-in_training,]

Budowa modelu

W wyniku przeprowadzonych badań model klasyfikacyjny został utworzony zgodnie z algorytmem Random Forest. Do budowy i opytymalizacji modelu wykorzystano powtórzoną ocenę krzyżową z podziałem na 2 części i 4 powtórzeniami.

regr_ctrl <- trainControl(
    # powtórzona ocena krzyżowa
    method = "repeatedcv",
    # liczba podziałów
    number = 2,
    # liczba powtórzeń
    repeats = 4)

rf_regr <- train(length ~ .,
             data = training,
             # Algorytm Random Forest
             method = "rf",
             trControl = regr_ctrl,
             ntree = 10,
             importance = TRUE)

Ostatecznie, w procesie uczenia jako optymalne zostały dobrane następujące wartości parametrów regresora.

## Random Forest 
## 
## 39438 samples
##    11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 4 times) 
## Summary of sample sizes: 19718, 19720, 19720, 19718, 19720, 19718, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared 
##    2    1.155766  0.5127557
##    6    1.123606  0.5393820
##   11    1.204624  0.4973097
## 
## RMSE was used to select the optimal model using  the smallest value.
## The final value used for the model was mtry = 6.

Ocena modelu

Oceny modelu dokonano na podstawie wcześniej przygotowanego zbioru testowego. Jako miary oceny przyjęto R2 oraz RMSE.

predicted <- predict(rf_regr, testing)
observed <- testing[, 'length']
  
SS_res <- sum( (observed - predicted)^2 )
SS_tot <- sum( (observed - mean(observed))^2 )
R_squared <- 1 - SS_res/SS_tot

RMSE <- sqrt( mean((observed - predicted)^2) / length(predicted) )
length(predicted)
## [1] 13144
miara wartość
R2 NA
RMSE 0.01

Ocena ważności atrybutów

Ocena ważności atrybutów ma za zadanie wskazać jakie atrybuty (czynniki) miały wpływ na to, że rozmiar śledzi zaczął w pewnym momencie maleć.

varImp(rf_regr)
## rf variable importance
## 
##        Overall
## X      100.000
## sst     51.647
## chel1   33.223
## cfin1   28.144
## chel2   23.050
## fbar    22.568
## totaln  16.474
## nao     14.276
## cfin2    6.206
## recr     5.410
## sal      0.000

Z powyższej tabeli wynika, że najistotniejszy jest atrybut sst reprezentujący temperaturę przy powierzchni wody. Nieco mniej ważny jest atrybut chel1 reprezentujący dostępność planktonu Calanus finmarchicus gat. 1. Warto zwrócić uwagę, że atrybut chel1 był silnie skorelowany (i z tego powodu odrzuconym z analizy, jako redundantny) atrybutem lcop1 (0.95). Zatem również dostępność planktonu w postaci widłonogów gat. 1 ma wpływ na zmniejszenie długości śledzi.

ToDo: Remove X parameter from analysis; correct data formating

Należy zauważyć, że istotność żadnego z atrybutów nie przekracza r, co może oznaczać, że także inne atrybuty nie uwzględnione w zbiorze danych mogły mieć równie istotny wpływ na zmniejszenie się długości śledzi. Przykładowe czynniki mogące mieć wpływ to zanieczyszczenie wody, czy dostępność światła.